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This paper considers a simple telephone exchange model which has an 
infinite number of trunks and in which the traffic depends on two parameters, 
the calling-rate and the mean holding-time. It is desired to estimate these 
parameters by observing the model continuously during a finite interval, 
and noting the calling-time and hang-up time of each call, insofar as these 
times fall within the interval. It is shown that the resulting information may, 
for the purpose of this estimate, be reduced without loss to four statistics. 
These statistics are the number of calls found at the start of observation, the 
number of calls arriving during observation, the number of calls terminated 
during observation, and the average number of calls existing during the 
interval of observation. The joint distribution of these sufficient statistics is 
determined, in principle, by deriving a generating function for it. From this 
generating function the means, variances, covariances, and correlation co- 
efficients are obtained. Various estimators for the parameters of the model 
are compared, and some of their distributions, means, and variances pre- 
sented. 

I THEORETICAL PROBLEMS AND METHODS OF TRAFFIC MEASUREMENT 

Four important kinds of theoretical problems arise in the measurement 
of telephone traffic. These are: (1) the choice of a mathematical model, 
containing parameters characteristic of the traffic, to serve as a descrip- 
tion; (2) the devising of efficient methods of estimating the parameters; 
(3) the determination of the anticipated accuracy of measurements; 
and (4) the assessment of actual accuracy, after measurements have 
been made. 

The present paper deals with aspects of the second and third kinds of 
problem, for the simplest and least realistic mathematical model of tele- 
phone traffic. Specifically, for this model, we treat the problems of (i) 
complete extraction of the information from a given observation period, 

939 



940 THE BELL SYSTEM TECHNICAL JOURNAL, JULY 1957 

without regard to costs of observation, and (ii) determination of the 
anticipated accuracy of certain methods of estimation which arise natu- 
rally from the discussion of complete extraction. 

The method by which we attack problems (i) and (ii) in this paper 
has three stages. First we choose a small number of significant properties 
of, or factors in, the physical system we are studying. Then we abstract 
these properties into a mathematical model of the physical system. Fi- 
nally, from the properties of the model, we derive results which may be 
interpreted as answers to the two problems treated. The advantage of 
this method is that we can use the precise, powerful apparatus of mathe- 
matics in studying the model; its limitation is that it yields results which 
are only as accurate as the model in describing reality. 

A method similar to the above forms the theoretical underpinning of 
telephone traffic engineering itself. To design equipment effectively, the 
traffic engineer needs a description of the traffic that is handled by central 
offices. He decides what properties of the entire system of telephone 
equipment and customers will be most useful to him in describing the 
traffic. He then designates certain parameters to serve as mathematically 
precise idealizations of these properties, and in terms of these parameters 
constructs a model of the traffic, upon which he bases much of his engi- 
neering. 

In choosing a mathematical model for a physical system, one is con- 
fronted with two generally opposed desiderata: fidelity to the system 
described, and mathematical simplicity. The model may involve impor- 
tant departures from physical reality; a model that is sufficiently amena- 
ble to mathematical analysis often results only after one has introduced 
admittedly" false assumptions, ignored certain effects and correlations, 
and generally oversimplified the system to be studied. However, the 
abstract model will be an exact and simple tool for analysis. 

We can construct a simple mathematical model for the operation of 
a telephone central office by leaving out of consideration many impor- 
tant facts about such systems, and by concentrating on factors most 
relevant to operation. Since we are interested in telephone traffic and 
in the availability of plant, it seems natural to require that a realistic- 
model take account of at least the following five significant factors: (1) 
the demand for telephone service; (2) the rate at which requests for 
service can be processed and connections established; (3) the lengths of 
conversations; (4) the supply of central office equipment; and (5) the 
manner in which the first four factors are interrelated. Unfortunately, 
the mathematical complexity of such a realistic; model precludes easy 
investigation. Therefore, the model used in this paper is based only on 
factors (1) and (3). 
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The demand for telephone traffic is usually made precise by describing 
a stochastic process which represents the way in which requests for tele- 
phone service occur in time. A realistic description will take account of 
the facts that, the demand is not constant, but has daily extremes, and 
that in small systems, the demand may be materially lessened when 
many conversations are in progress. Since taking account of the first fact 
leads to a more complicated model in which our investigations are more 
difficult, we ignore it, with the proviso that the results we derive are 
only applicable to systems and observations for which the demand is 
nearly constant. The second kind of variation in demand becomes insig- 
nificant as the number of subscribers increases and the traffic remains 
constant. Hence, we further confine the applicability of our results to 
systems with large numbers of subscribers, and we assume that the de- 
mand does not depend on the number of conversations in existence. 

With these assumptions, a mathematically convenient description of 
the demand is specified by the condition that the time-intervals between 
requests for service have lengths which are mutually independent posi- 
tive random variables, with a negative exponential distribution. 

A telephone central office contains two kinds of equipment: control 
circuits which establish a desired connection, and talking paths over 
which a conversation takes place. The time that a request for service 
occupies a unit of equipment, be the unit a control circuit or a talking 
path, is called the holding-time of the unit. A request for service affects 
the availability of both kinds of equipment but, except for special cases, 
the holding-times of talking paths are usually much longer than the 
holding-times of control units such as markers, connectors, or registers. 
In view of this disparity, we assume that the only holding-times of con- 
sequence are the lengths of conversations; i.e., the holding-times of 
talking paths. We assume also that these lengths are mutually inde- 
pendent positive random variables, with a negative exponential distribu- 
tion. 

For the simplest mathematical model of telephone traffic, we may 
consider the arrangement of switches and transmission lines which con- 
stitutes a talking path in the physical office to be replaced by an abstract 
unit called a "trunk". A trunk is then an abstraction of the equipment 
made unavailable by one conversation, and we may measure the supply 
of talking paths in the office by the number of trunks in a model. The 
word "trunk" is also used to mean a transmission line linking two central 
offices, but as long as we have explained our use of the word there need 
be no confusion. Often the number of transmission lines leading out of 
an office is a major limitation on its capacity to carry conversations, 
and in this case the two uses of the word "trunk" are very similar. Un- 
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fortunately, we do not take advantage of this similarity, since we make 
the mathematically convenient but wholly unrealistic assumption that 
the number of trunks in the model is infinite. 

The model we investigate thus depends on only two of the factors 
previously listed as essential to a realistic model: namely, (1) the demand 
for service, and (3) the lengths of conversations. In view of the simplicity 
and inaccuracy of this model, the question arises whether much is gained 
from a detailed analysis. Such scrutiny may indeed reveal little that is 
of great practical value to traffic engineers. It is important methodologi- 
cally, however, to have a detailed treatment of at least one approximate 
case. We undertake this detailed treatment largely for the insight that 
it may give into methods which could be useful in dealing with more 
complex and more accurate models. 

Once a designer has chosen a model and has specified the parameters 
he would like to have measured, it is up to the statistician to invent effi- 
cient means of measurement, by choosing, for each parameter, some 
function of possible observations to serve as an estimate of that parame- 
ter. One measure of efficiency that is of mostly theoretical interest is the 
observation time required to achieve a given degree of anticipated ac- 
curacy; the most realistic measure of efficiency is in terms of dollars and 
man-hours. It may often be more efficient, in the sense of the latter 
measure, to spread observation over enough more time to compensate 
for the inability of an intrinsically cheaper method of measurement to 
extract all of the information present in a fixed time of observation. For 
example, periodic scanning of switches in a telephone exchange is usually 
less costly than continuous observation. As a result, telephone traffic 
measurement is usually carried out by averaging sequences of instan- 
taneous periodic observations of the number of calls present, rather than 
by continuous time averaging, although it can be shown that continuous 
observation is more efficient at extracting information. Thus statistical 
efficiency, which may be expensive in terms of measuring equipment, 
can be exchanged for observation time, which may be less costly. This 
exchange brings about a reduction in cost without impairing accuracy. 

Our concern in this paper is with the less practical problems of com- 
plete extraction, and of the anticipated accuracy of estimation methods 
based on complete extraction. Let us consider how our mathematical 
model can shed light on these problems. A mathematical model may or 
may not be a faithful description of the behavior of real telephone sys- 
tems. Nevertheless random numbers, with or without modern computing 
machines, enable one to make experiments and observations on physical 
situations which approximate, arbitrarily closely, any mathematical 
model. Thus we can speak meaningfully of events in the model, and of 
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making measurements and observations on the model. The mathematical 
model elucidates our problems in the following ways: (1) it enables us 
to state precisely what information is provided by observation; (2) it 
enables us to explain what we mean by complete extraction of informa- 
tion; and (3) it enables us to derive results about the anticipated ac- 
curacy of measurements in the model. These results will have approxi- 
mately true analogues in physical situations to which the model is 
applicable. 

The calls existing during the observation interval (0, T) fall into four 
categories: (i) those which exist at 0, and terminate before T; (ii) those 
which fall entirely within (0, T) ; (iii) those which exist at and last 
beyond T; and (iv) those which begin within (0, T) and last beyond T. 
For calls of category (i), we assume that we observe the hang-up time 
of each call; for category (ii), we observe the matching calling-time and 
hang-up time of each conversation; for category (iii), we observe simply 
the number of such calls; and for category (iv), we observe the calling- 
times. Table I summarizes the kinds of calls and the information ob- 
served about each. 

What we mean by the complete extraction of information is made 
precise by the statistical concept of sufficiency. By a statistic we mean 
any function of the observations, and by an estimator we mean 
a statistic which has been chosen to serve as an estimate of a particular 
parameter. Roughly and generally, a set S of statistics is sufficient for a 
set P of parameters when S contains all the information in the original 
data that was relevant to parameters in P. If S is sufficient for P, there 
is a set E of estimators for parameters in P, such that the estimators in 
E depend only on statistics from S, and such that an estimator from E 
does at least as well as any other estimator we might choose for the same 
parameter. Thus we incur no loss in reducing the original data (of speci- 
fied form) to the set S of statistics. It remains to state what it means for 
S to contain all the relevant information. We do this in terms of our 
model. 

The mathematical model we are adopting contains two distribution 





Table I — Information 


Observed 


Types of Calls 


Start in (0, T) 


Start before 


End in (0, T) 


(ii) , matching calling-times and 
hang-up times known, num- 
ber of calls known 


(i) , hang-up times known, num- 
ber of calls known 


End after T 


(iv), calling-times known, num- 
ber of calls known 


(iii), number of calls known 
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functions, that of the intervals between demands for service, and that 
of the lengths of conversations. We have supposed that these distribu- 
tions are both of negative exponential type, each depending on a single 
parameter. Thus we know the functional form of each distribution, and 
each such form has one unknown constant in it. Since the mathematical 
structure of the model is fully specified except for the values of the two 
unknown constants, we can assign a likelihood or a probability density 
to any sequence 2 of events in the model during the interval (0, T). 
This likelihood will depend on the parameters, on 2, and on the number 
of calls in existence at the start of the interval. If the likelihood L(2) 
can be factored into the form L = FH, where F depends on the param- 
eters and on statistics from the set S only, and H is independent of the 
parameters, then the set S of statistics may be said to summarize all the 
information (in a sequence 2) relevant to the parameters. If L can be 
so factored, then S is sufficient for the estimation of the parameters. 

The mathematical model to be used in this paper is described and 
discussed in Sections II and III, respectively. Section IV contains a 
summary of notations and abbreviations which have been used to sim- 
plify formulas. 

In Appendix A we show that the original data we have allowed our- 
selves can be replaced by four statistics, which are sufficient for estima- 
tion. In Appendix B and Sections V-VIII we discuss various estimators 
(for parameters of the model) based on these four statistics. To determine 
the anticipated accuracy of these methods of measurement, we consider 
the statistics themselves as random variables whose distributions are 
to be deduced from the structure of the model. 

A primary task is the determination of the joint distribution of the 
sufficient statistics. In view of the sufficiency, this joint distribution tells 
us, in principle, just what it is possible to learn from a sample of length 
T in this simple model. By analyzing this distribution we can derive 
results about the anticipated accuracy of measurements in the model. 

The joint distribution of the sufficient statistics is obtainable in prin- 
ciple from a generating function computed in Appendix C, using methods 
exemplified in Section X. This generating function is the basic result of 
this paper. The implications of this result are summarized in Section 
IX, which quotes the generating function itself, and presents some 
statistical properties of the sufficient statistics in the form of four tables: 
(i) a table of generating functions obtainable from the basic one; (ii) a 
table of mean values; (iii) a table of variances and covariances; and 
(iv), a table of squared correlation coefficients. (The coefficients are 
all non-negative.) 
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II DESCRIPTION OF THE MATHEMATICAL MODEL 

Throughout the rest of the paper we follow a .simplified form of the 
notational conventions of J. Riordan's paper 11 wherever possible. A sum- 
mary of notations is given in Section IV. The model we study has the 
following properties: 

(i) Demands for service arise individually and collectively at random 
at the rate of a calls per second. Thus the chance of one or more demands 
in a small time-interval At is 

aAt + o(At), 

where o(At) denotes a quantity of order smaller than At. The chance of 
more than one demand in At is of order smaller than At. It can be shown 
(Feller, 2 p. 364 et seq.) that this description of the demand is equivalent 
to saying that the intervals between successive demands for service are 
all independent, with the negative exponential distribution 

1 — e . 

This again is equivalent to saying that the call arrivals form a Poisson 
process; 2 i.e., that for any time interval, t, the probability that exactly n 
demands are registered in t is 

c- at (at) n 
n ! 

Thus the number of demands in / has a Poisson distribution with mean 
at. 

(ii) The holding-times of distinct conversations are independent vari- 
ates having the negative exponential distribution 

i - r", 

where y is the reciprocal of the mean holding-time h. This description of 
the holding-time distribution is the same as saying that the probability 
that a conversation, which is in progress, ends during a small time- 
interval At is 

yAt + o(At), 

without regard to the length of time that the conversation has lasted 
Feller, p. 375). 

(iii) The model contains an infinite number of trunks. Thus, at no time 
will there be insufficient central office equipment to handle a demand 
for service, and no provision need be made for dealing with demands that 
cannot be satisfied. 
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The original work on this particular model for telephone traffic is in 
Palm, 9 and Palm's results have been reported by Feller 3 and Jensen. 
The results have been extended heuristically to arbitrary absolutely con- 
tinuous holding-time distributions by Riordan, 11 following some ideas of 
Newland 8 suggested by S. 0. Rice. 

Let Pn(t) be the probability that there are j trunks busy at t if there 
were i busy at 0. And let P,-(f, x) be the generating function of these 
probabilities, defined by 

Pi(t,x) = Es'PtfCfl. 

;=0 

Then Palm 9 has shown (pp. 56 et seq.) that 

Pi(t, x) - [1 + Or - 1) e- yt Y exp {(x - l)ah (1 - e"")}. 

This is formula (12) of Riordan 11 with his g replaced by e~ yt . It can be 
verified that the random variable N(t) is Markovian; the limit of P,(/, x) 
as t — * *> is 

exp {(.r - I) ah), 

so that the equilibrium distribution of the number of trunks in use is a 
Poisson distribution with mean b = ah. The shifted random variable 
[N(t) — b] then has mean zero, and covariance function be"' . 

For additional work on this model the reader is referred to F. W. 
Rabe, 10 and to H. Stormer. 12 

Ill DISCUSSION OF THE MODEL 

Let us envisage the operation of the model we have described by con- 
sidering the random variable N(t) equal to the number of trunks busy 
at time t. As a random function of time, N(t) jumps up one unit step each 
time a demand for service occurs, and it jumps clown one unit step each 
time a conversation ends. If N(0 reaches zero, it stays there until there 
is another demand for service. If N(t) = n, the probability that a con- 
versation ends in the next small time-interval At is 

nyAt 4- o(M), 

because the n conversations are mutually independent. A graph of a 
sample of N{t) is shown in Fig. 1. 

The model we described departs from reality in several important 
ways, which it is well to discuss. First, the assumption that the number 
of trunks is infinite is not realistic, and is justified only by the mathe- 
matical complication which results when we assume the number of trunks 
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TIME, t 



Fig. 1 — A graph of N(t). 

to be finite. It can also be argued that unlimited office capacity is ap- 
proached by offices with adequate facilities and low calling rates, and 
therefore, in some practical cases at least, the model is not flagrantly 
inaccurate. 

Second, the choice of a constant calling rate for the model ignores the 
fact that in most offices the calling rate is periodic. Thus, the applica- 
bility of our results to offices whose calling rates undergo drastic changes 
in time is restricted to intervals during which the normally variable 
calling rate is nearly constant. Finally, although the assumption of a 
negative exponential distribution of holding-time affords the model great 
mathematical convenience, it is doubtful whether in a realistic model 
the most likely holding-time would have length zero, as it does in the 
present one. 

IV SUMMARY OF NOTATIONS 

a = Poisson calling rate 

// = mean holding-time 

7 = /r 1 = hang-up rate per talking subscriber 

b = ah = average number of busy trunks 
N(f) = number of trunks in use at t 
(0, T) = interval of observation 

n = N(0) = number of trunks in use at the start of observation 
A = number of calls arriving in (0, T) 
H = number of hang-ups in (0, T) 
K = A+H 
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Z = ( N(t) 
Jq 



M = Z/T = average of N(t) over (0, T) 

\p„] = the (discrete) probability distribution of n, the number of 
trunks found busy at the start of observation 

An estimator for a parameter is denoted by adding a cap ( A ) and a 
subscript. The subscripts differentiate among various estimators for the 
same parameter. We use d e = A/T, y c = H/Z, di = K/2T, 
71 = K/2Z, and 72 = A/Z. 

Also, it is convenient to use the following abbreviations: r for yT, and 
C for (1 — e~ r )/'r, where r is the dimensionless ratio of observation-time 
to mean holding-time. The symbol E is used throughout to mean mathe- 
matical expectation. 

V THE AVERAGE TRAFFIC 

We have adopted a model which depends on two parameters, the 
calling rate a, and the mean holding-time h, or its reciprocal 7. Before 
searching for a set of statistics that is sufficient for the estimation of 
these parameters, let us consider the product ah = b. This product is 
important because, as we saw in Section II, the equilibrium distribution 
of the number of trunks in use depends only on b, and not on a and h 
individually. Indeed, the equilibrium probability that n trunks are busy 
is 

— bin 

e b 



n 



and the average number of busy trunks in equilibrium is just b. 
The average number of trunks busy during a time interval T is 



M - i / N{t) dl; 



i.e., the integral of the random function N(t) over the interval T, divided 
by T. This suggests that for large time intervals T, M will come close 
to the value of b, and can be used as an estimator of 6. Since M is a ran- 
dom variable, the question arises, what are the statistical properties of 
Ml This question has been considered in the literature, the principal 
references being to F. W. Rabe 10 and to J. Riordan. 11 Riordan's paper is 
a determination of the first four semi-invariants of the distribution of M 
during a period of statistical equilibrium, but without restriction on the 
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assumed frequency distribution of holding-time. It follows from Gior- 
dan's results that M converges to b in the mean, which is to say that 

lim E { \M - b | 2 } = 0. 

T-x 

It also follows that M is an unbiased estimator of 6; i.e., that E\M\ = b, 
and that M is a consistent estimator of b, which means that 



lim pr{\ M — b \ > e] = 

T-.00 



for each e > 0. 



VI MAXIMUM CONDITIONAL LIKELIHOOD ESTIMATORS 

As shown in Appendix A, the likelihood L c of an observed sequence, 
conditional on N(0), is defined by 

In L e = A In a + H In y - yZ - aT. 

According to the method of maximum likelihood, we should select, as 
estimators of a and y respectively, quantities d ( . and 7,. which maximize 
the likelihood L, . Now a maximum of L c is also one of In L , and vice 
versa. Therefore a c and 7,. are determined as roots of the following two 
equations, called the likelihood equations: 

AlnL t = 0; f\nL c = 0. 
da 07 

The solutions to the likelihood equations are 

A A _H 

a <- ~ Ipy Tc — y • 

These are the maximum conditional likelihood estimators of a and 7. 
The estimator d c is the number of requests for service in T divided by T; 
this is intuitively satisfactory, since d e estimates a calling rate. 

Since maximum likelihood estimators of functions of parameters are 
generally the same functions of maximum likelihood estimators of the 
parameters, we see that AZ/HT is a maximum likelihood estimator of b. 

VII PRACTICAL ESTIMATORS SUGGESTED BY MAXIMIZING THE LIKELIHOOD 

L, DEFINED IN APPENDIX A 

We obtain as likelihood equations 

A in L - 0, ~ In L = 0. 

do ay 
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These may be written as 



and 



n + A 



7 = 1 

z + n 

7 



The first of these shows the estimated calling rate as a pooled combination 
of the conditional estimate A/T, considered in the last section, and an 
estimate n/h based on the initial state. This latter estimate has the form 

calls in progress 



mean holding time' 



and so is intuitively reasonable, since b/h = a. The second equation 
exhibits our estimate of 7 as a pooled combination of the conditional 
estimate H/Z and the ratio a/n. This ratio is acceptable as an estimate 
of 7, since a/b = 7 and b — E{n] is the average value of n. 

If we substitute, in the right-hand sides of these equations, the condi- 
tional estimators A/T, H/Z, and Z/H for o, 7, and h, respectively, we 
obtain simple, intuitive estimators which include the influence of the 
initial state n, and show how it decreases with increasing T. Thus 



n + A 


estimates a, 




H + TH 

7 ,nZ 
Z + -jj 


estimates 7. 



VIII OTHER ESTIMATORS 



Additional estimators may be arrived at by intuitive considerations, 
or by modifying certain maximum likelihood estimators. Some estimators 
so obtained are important because they use more of the information 
available in an observation than do the conditional estimators d c and 
7 C , without being so complicated functionally that we cannot easily 
study their statistical properties. 
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It seems reasonable, and can be shown rigorously (Appendix C), that 
for an interval (0, T) of statistical equilibrium, the distribution of .4 
and that of H are the same. Thus we can argue that, for long time inter- 
vals, A and H will not be very different. This suggests using 

A _ A + H = K_ 

1 ~ 2T 2T 

as an estimator of a. This estimator does not involve 7, and it uses not 
only information given by A , but also information supplied by arrivals 
occurring possibly before the start of observation. 

Similarly, since b = a/7, and M is a consistent and unbiased estimator 
of b, we may use 

" = *L -= L 

71 2Z J, 

to estimate 7, and its reciprocal to estimate h. Finally, since for long (0,T) 
we have A ~ H, we may try 

A 1 

72 = 7 = I 

as an estimator of 7, and its reciprocal as an estimator of h. 

IX THE JOINT DISTRIBUTION OF THE SUFFICIENT STATISTICS 

The basic result of this paper is a formula for the generating function 

E{z n x MT) w A u H e- iz ) (9.1) 

for the joint distribution of the random variables n, N(T), A, H, and Z. 
This formula is derived in Appendix C, by methods illustrated in Section 
X. For an initial n distribution \p n ) , the generating function is 

"(M: + yx - 7 ")<T (m)r + yuT 



;i>0 



f -' 7 ^ (9.2) 



' exp \ (r + 7) 2 " + 7T7 " aT 

It is proved in Appendix A that the set of statistics {n, A, H, Z\ is 
sufficient for estimation on the basis of the information assumed, which 
was described in Section I. Thus the generating function (9.2) specifies, 
at least in principle, what can be discovered about the process from an 
observation interval (0, T), for which N(0) has the distribution \p n }. 
All the results summarized in this section are consequences of (9.2). 
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Table II 



In £(A') 



1. e-f z 

2. e-f^ 

3. >J K 

4. e"f°. 

5. ijKe-i* 



b -f T + 

b\-i + 



f 2 T r 2 (l - e-«+t )T )' 



r + 7 (f + 7) 2 

r 2 _ f 2 (l - e-<f+ r >) ' 
f + r (f + r) J 

2orC(y - 1) + ar(l - Ofo 1 - 1) 
2aTC(e-f /2r - 1) + oT(l - C)(e-f'r - 1) 



'(-Ff-J 



e -f( + r) _ 1] _ r f 1 - 



*)] 



By substitution, and by either letting the appropriate power series 
variables — > 1, or letting f — > 0, or both, we can obtain from (9.2) the 
generating function of any combination of linear functions of the basic 
random variables n, N(T), A, H, and Z. Some of the generating func- 
tions thereby obtained are listed in Table II, in which the entries all 
refer to an interval (0, T) of equilibrium. 

Since, for equilibrium (0, T), the generating functions are all exponen- 
tials, it has been convenient to make Table II a table of logarithms of 
expectations, with random variables X on the left, and functions In 
E{X) on the right. C as a function of r is plotted in Fig. 2. 

Entry 1 of Table II is actually the cumulant generating function of 
Z for equilibrium (0, T) ; similarly, Entry 2 is that of M, and depends 
only on the average traffic b and the ratio r. The forni of the general 
cumulant of M is 



I " Jo 



dx. 



This result coincides with a special case (exponential holding-time) of 
a conjecture of Riordan." This conjecture was first established (for a 
general holding-time distribution) in unpublished work of S. P. Lloyd. 
The cumulant generating function permits investigation of asymptotic 
properties. We prove in Section X that the standardized variable 

v = (yT/Zb) m (M - b) 

= (r/2b) in (M - b) 

is asymptotically normally distributed with mean and variance 1. 
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i.o 



'C 



0.4 



0.2 



c-qP 



12 



Fig. 2 — C as a function of r. 

From Entry 3 of Table II it can be seen that K is distributed as 
2w + v, where u and u follow independent Poisson distributions with the 
respective parameters aT{\ — C) and 2aTC. The probability that K = 
n for an interval of equilibrium is 



r„ = exp [aT{C - 1)} £ 



(2arC)"~ 2i (of - aW 
(n-2j)! j| 



where the sum is over fs for which 0^2^'^ n. 

The estimator <2i for a is equal to K/2T, and has mean and variance 
given by 

E\di} = o, 



var [di] - ^ (2 - C). 

The distribution of «i is given by 

pr{«i g x] = E r »> 

the summation being over n ^ 2T.r. 

From (9.2) one can obtain, by substitution of the stationary n distribu- 
tion for {p n }, and subsequent differentiation, the means, variances, 
covariances, and correlation coefficients of the sufficient statistics, for 



<)51: 
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Table III — E{X, Y] 





1 


n 


A 


H 


K 


Z 


1 


1 


I, 


aT 


aT 


2aT 


bT 


n 




6(1 + 6) 


baT 


aT(,C + 6) 


aT(C + 26) 


bT(C + b) 


A 






aT(l + aT) 


aT(l - C+aT) 


aZ*(2 -C + aT) 


bT(\ -C + aT) 


If 








aT{\ + aT) 


aT(2 -C + aT) 


bT{\ -C + aT) 


K 










2aT(2 - C + 2aT) 


26T(1 -C + aT) 


Z 












bTh{2 -2C + aT) 



Table IV — cov{X,Y} 





n 


A 


H 


K 


Z 


n 
A 
H 
K 
Z 


b 




aT 


aTC 

aTi\ - C) 
aT 


aTC 

aT(2 - C) 
aT(2 - C) 
2aT{2 - C) 


bTC 

bT(l - C) 
bT{\ - C) 
2bT(l - C) 
2bTh{\ - C) 



equilibrium intervals (0, T). It has been convenient to display these in 
three triangular arrays, the first consisting of expectations of products, 
the second comprising the variances and covariances, and the third 
exhibiting, for simplicity, the squared correlation coefficients, since the 
correlation coefficients are never negative for these random variables. 

In Table III, the entry with coordinates (X, Y) \sE\XY\ for equilib- 
rium (0, T). All three tables are expressed in terms of a, b, T, h, r, and 
C, the last of which is plotted in Fig. 2. 

The variances and covariances of the sufficient statistics are listed in 
Table IV; the entries are of the form: 

cov {X, Y\ - E{XY) - E\X\E{Y\. 

Table V, finally, lists the squared correlation coefficients; i.e., the 
quantities 

cov 2 \X,Y) 



p\X, Y) = 



var {X\ var \Y)' 



For any time interval (0, T), A has a Poisson distribution with param- 
eter aT, so that Td c does also. Therefore the distribution of d c is given 

by 



pr \d c ^ x) = X) 



r (aT)" 



nl 



where the summation is over n ^ xT. Evidently 



and 
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E{d c } = a, 



var \a c \ = ^, 



V~>» 



.so that d c is an unbiased and consistent estimator of a. We now compare 
the variances of estimators d r and d\. From Table IV we have 



var 



l«i) = Jfl - g ) < % = var (*■!» 



so that di is a better estimator of a for any T > 0, in the sense that its 
variance is less. 



X THE DISTRIBUTIONS of Z AND M 



Since we have defined 



Jo 



iV(0 df, 



we can regard Z as the result of growth whose rate is given by the ran- 
dom step-function N(t); when N(t) = n, Z is growing at rate n. An idea 
similar to this is used by Kosten, Manning, and Garwood", and by Kos- 
ten alone. 5 Now the Z(T) process by itself is not Markovian, but it can 
be seen that the two-dimensional variable {N(l), Z(t)\ itself is Marko- 
vian. Let F n (z, I) be the probability that N(t) = n and Z(t) ^ z. Since 
the two-dimensional process is Markovian, we can derive infinitesimal 
relations for F n (z, t) by considering the possible changes in the system 
during a small interval of time Af. 

Table V — p 2 (X ,Y) 





u 


,1 


H 


A' 


Z 




1 





1 - e~ r 


rC 2 
2 - C 


rC* 




2(1 - C) 


A 




1 


1 - C 


2 - C 
2 


1 - C 
2 


H 






1 


2 - C 


1 - C 








2 


2 


K 








1 


1 - c 

2 - C 


Z 










1 
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If N(t) = n, then the probability is [1 - ynAt - a At - o(M)] that 
there is neither a request for service nor a hang-up during At following 
/, and that Z{t + At) = Z(t) + ?iAt. Therefore the conditional proba- 
bility that N(t + At) = n and Z(t + At) ^ z, given that no changes 
occurred in At, is 

F n (z - nAt, t). 

For N(t) = {n + 1), the probability is y(n + l)At + o(At) that one 
conversation will end during At following t. The increment to Z(t) 
during At will depend on the length x of the interval from / to the point 
within At at which the conversation ended. The increment has magni- 
tude (n + l).r + ?i(At - x) = x + ft At, as can be verified from Fig. 3, 
in which the shaded area is the increment. Since x is distributed uniformly 
between and At, the increment x + nAt is distributed uniformly be- 
tween nA2 and (ft + 1)A£. Therefore the conditional probability that 
N(t + At) = ft and Z(t + At) ^ z, given that one conversation ended 
in At, is 

.(n+l)A( 

F n+1 (z — ft, t) dn. 



At JnAJ 

By a similar argument it can be shoAvn that the probability that one 
request for service arrives in At is aAt + o(At), and that the conditional 
probability that N{t + At) = n and Z(t + At) ^ z, given that one 
request arrived during At, is 

I r" At 

— / F„_i(z - ft, t) du. 

At J(n-1)A( 

Define F„(z, t) to be identically for negative ft. Adding up the probabil- 



n+2 - 
n+ i 
n 
n-i 
n-2 




Fig. 3 — Increment to Z in At. 
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ities of mutually exclusive events, we obtain the following infinitesimal 
relations for F„(z, /): 

*(n4-l)Al 

F„(z, t + At) - y(n + 1) F n+1 (z - u, t) du 

+ a / F„_,(z - u, t) du + F„(z - nAt, t) 

■/(»— 1)A« 

•fl — A<(7« + o)] + o(A/), for any n. 

Expanding the penultimate term of the right side in powers of nAt, 
and the left side in powers of At, we divide by At, and take the limit as 
At approaches 0. Now 

lim — / F n+1 (z - u, t) du = F„ +l (z, t). 

A*-0 At J nM 

Thus, omitting functional dependence on z and / for convenience, we 
reach the following partial differential equations for F„(z, t): 

ir.-^ + i)r^ + a r^-n±F. (101) 

— [yn + a]F„ , for any n. 

Since Z(0) = 0, Ave impose the following boundary conditions: 

F„(0, 0=0 for n > and t > 0, 

F n (z, 0) = p n for z = 0, (10.2) 

F„(z, 0) = for z < 0, 

where the sequence \p„] forms an arbitrary N(0) distribution that is 
zero for negative n. 

To transform the equations, we introduce the Laplace-Stieltjes in- 
tegrals 

*.(f, t) = \ e'" dF n (z, t), / ^ 0, Re (f) > 0, 
Jo- 
in which the Stieltjes integration is understood always to be on the 
variable z. We note that 



and that 



(" e-"F n (z,t)dz =^„(f,/), 
Jo- J 

Vn (t, t) - Fn(0, i) + [ a"* | F„ (z, dz. 

Jo oz 



958 THE BELL SYSTEM TECHNICAL JOURNAL, JULY 1957 

Applying now the Laplace-Stieltjes transformation to (10.1), we obtain 
J£ = y( n + 1V„ +1 + cup„-i - nfa, + n£F n (Q, t) ^ Q ^ 

- [7*1 + a]<p„ , 

in which we have left out functional dependence on f and / where it is 
unnecessary. By the boundary conditions (10.2), n{F n {0, t) = for 
n ^ and / > 0; in (10.3) we may therefore omit this term in the region 
/ > 0. Let <p be defined by 

;i=0 

The series is absolutely convergent for \x\ < 1 , since 

|«>„(r,*) I ^ 1, for all n. 
The following partial differential equation for <p is obtained from 
(10.3): 

% + [f.r + tCt - 1)1 d / = a(x - IV. (10.4) 

If we integrate out the information about Z by letting £ approach in this 
equation, we obtain the equation derived by Palm (loc. cit.) for the gene- 
rating function of N(t). Therefore our equation has a solution of the 
same form as Palm's. For the boundary conditions (10.2), this solution is 

exp { — - — ; — ^ — Ifa + TU - 1)1 - 



(f + - )! f + 7 ' (loo 

,=. fit * + 7(x - 1)l<T" +r " + 7 T 

•S p ". FF^ J' 

Actually <p contains more information than we want since it yields the 
joint distribution of N and Z. We may integrate out the former variable 
by letting x approach 1 in 10.5. Then, 

£{exp (-rZ)I = exp ( U (r + 7)2 y^ 



■£*P 



-<r+7)r 






is the Laplace transform of the distribution of Z for an arbitrary N(0) 
distribution \p„\. This result is not restricted to an interval (0, T) 



STATISTICS FOR A SIMPLE TELEPHONE EXCHANGE MODEL 959 

of statistical equilibrium; however, if the sequence \p„\ does form the 
stationary distribution discussed in Section II, then 

X»„ = cxp{b(x- 1)), (10.6) 

and 

+ = exp 1 (r + 7) 2 ~ V^-y) (107) 

is the Laplace transform of the distribution of Z for an interval (0, T) 
of statistical equilibrium. 
The Laplace transform is a moment generating function expressible as 

, = y> (-f)"w n 
* fe nl ' 

where m, is the rc" 1 ordinary moment of Z. Differentiation of 10.7 then 
gives a recurrence relation for the moments upon equating powers of 

(-r). Thus, 

= ^-{2ar(l - e-« +y)T ) + (f + 7)[7«7' + &^- (f+T,T ]), 
and 
7 3 w„ + i - 3y 2 nm n + 3-yn(n - l)wi„-i - n(n — l)(w - 2)m„_ 2 

= rfftn. - (2a + ayT)nmn-i + 2ane~ TT (»i + T)"~ l + n (10.8) 
•(» -l)aTe- yT (m + r)" -2 - »(n - l)(w - 2)bTe~' lT {m + T)"-\ 
where (ra + T) n is the usual symbolic abbreviation of 



§(?)*-- 



From the recurrence (10.8) it is easily verified that 
rni = bT, 

m = (bT) 2 + — [i - c], 

7 
from which it follows that the variance of Z is 

9hT 

var [Z\ = — [1 - C]. 

7 



9G0 
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Since \f/ is the Laplace-Stieltjes transform of the distribution of Z 
over an interval of equilibrium, In \p is the cumulant generating function, 
and has the following simple form : 



In ^ = 6 



%t„-U+y'> T 



f(e 



- 1) 



(f + t) 2 



..r- 



?t + 






ytT 1 

f»(l - e -tf+T)»-) 

(f + 7) 2 



]■ 



(10.9) 



M is a linear function of Z, so we may obtain the cumulant generating 
function of M in accordance with CrameV (p. 187). This function is 



rr*-t\~l 



b\ -f + 



r + f 



(r + f) 2 



(10.10) 



and depends only on 6 and r. 

The mean and variance of M for an interval of equilibrium are respec- 
tively given by 



E[M] = b, 
var {M } = - [1 - C], with C = 



1 - e" 



results which were first proved in Riordan." A normal distribution having 
the mean and variance of M has the cumulant generating function 



b r_ r + £ + «vl- o 

r r 2 



]■ 



(10.11) 



which is to be compared to (10.10). Since var {M } goes to as T ap- 
proaches oo , we may expect that a suitably normalized version of Z will 
be asymptotically normally distributed as T approaches °° . The cumu- 
lant generating function of the normalized variable (2bhT)~ ll2 (Z — bT) 

is 



- r>- 1 



•&) 



1/2 



+ 2 



exp(-rg) '->•} 



which approaches f 2 / 2 as T —*■ °o . It follows that the normalized variable 
is asymptotically normal with mean and variance 1, and that 
(r/2b) h (M — b) is also asymptotically normal (0, 1). 
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Appendix A 

PROOF THAT {n, A,H, Z\ IS SUFFICIENT. 

We observe the system during the interval (0, T), and gather the in- 
formation specified in Section I, and summarized in Table I. From this 
information we can extract four sets of numbers, described as follows: 
S a the set of complete observed inter-arrival times, not counting 

the interval from the last arrival until T 
Sh the set of complete observed holding times 
Si the set of hang-up times for calls of category (i) 
Si the set of calling-times for calls of category (iv) 
In addition, our data enable us to determine the following numbers: 
n the number N{0) of calls found at the start of observation 
k the number of calls of category (iii); i.e., of calls which last through- 
out the interval (0, T) 
x the length of the time-interval between the last observed arrival 

and T 
In view of the negative exponential distributions which have been 
assumed for the inter-arrival times and the holding-times, and in view 
of the assumptions of independence, we can write the likelihood of an 
observed sequence of events as 

L—kyT—az TT —au TT —f- TT —yw TT — y(T— y) 

= e Vn • 11 ae ■ H ye y ■ H e H e w , 

uts a -<' s '/i u't.Si ytSt 

so that 

In L = —ykT — ax + In p n + A In a — ^ au 

utS a 

+ H In 7 - S 72 - St^-S y(T — y) 

zes h wcSi ytSi 

It is easily seen that the summations and the two initial terms can be 
combined into a single term, so that we obtain 

In L = In p n + A\na + H\ny -yZ - aT. 

This shows that L depends only on the statistics n, A, H, and Z; it 
follows that the information we have assumed can be replaced by the 
set of statistics \n, A, H, Z), and that these are sufficient for estimation 
based on that information. 

The likelihood is sometimes defined without reference to the initial 
state, by leaving the factor p„ out of the expression for L. Strictly speak- 
ing, this omission defines the conditional likelihood for the observed 
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sequence, conditional on starting at n. We use the notation: 

Pn 

A definition of likelihood as L c has been used by Moran. 7 Clearly 
In L c = A In a + H In y - yZ - aT. 

Appendix B 

UNCONDITIONAL MAXIMUM LIKELIHOOD ESTIMATES 

The definition of likelihood as L leads to complicated results which 
are of theoretical rather than practical interest. For this reason these 
results have been relegated to an appendix. 

The results of setting d/dy In L and d/da In L equal to zero lead, re- 
spectively, to the likelihood equations 

a — y(n — H) — y Z = 0, 

yn — a + yA — ayT = 0. 

Considered as a system of equations for y and a, this pair has the non- 
negative roots 

H - n - M + {(H - n - M) 2 + 4MK} 1/2 



7 = 



a = y - yM. 



2Z 



These are the unconditional maximum likelihood estimators for y and a. 
Although d c depended only on A and T, and y c only on H and Z, the 
unconditional estimators depend on all of n, A , H, Z, and T. We may 
obtain a maximum unconditional likelihood estimator for b as well, 
either by considering L to be a function of b and y, or from general 
properties of maximum likelihood estimators. Since b = a/y, we expect 
that b = a/y, as can be verified by an argument similar to that used 
above for & and y . 

The estimators a, b, and y obtained in this Appendix may turn out to 
be useful in practice, but their complicated dependence on the sufficient 
statistics n, A, H, and Z makes a study of their statistical properties 
difficult. As a first step along such a study, we have derived the gen- 
erating function of the joint distribution of the sufficient statistics in 
Appendix C. The greater simplicity of the conditional estimators of 
Section VI makes it possible to study their statistical properties. This 
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fact gives them a practical ascendancy over the unconditional estimators, 
even though the latter may be more efficient statistically by dint of 
using all the information available in an observation. 

Appendix C 

THE JOINT DISTRIBUTION OF N(t), 71, A, H, AND Z 

By methods already used in Section X one can obtain a gen- 
erating function for the joint distribution of all the random variables n, 
N(t), A, H, and Z. Let 

i Rf| S(t) A B —tZ\ 

& = E \x w u e \. 

Then <i> satisfies the differential equation 

d<i> d<p 

-TT + fax + 73! - yu] — = a(wx - 1)$, 

at ox 

whose solution has the form 

* = R{\$x +7* - yu]e- lt+y)i ) 

(aw\$x + yx - yu][l - e~ ((+y)t \ . ayvmt \ 

• exp v w+y? + tt-j ~ at r 

where the function R is determined by the initial distribution {p n \ 
through the relation 



8U1 _ v „ rt+2-T. 

n^O Lf + 7 J 



From these results it follows that the generating function 

<j\z x w u e ' } 
is given by 

v ■ f (<;x+yx-yu)e-« +y)T +yu \ n 
»i?o \ f + 7 / 

( aw(lx +yx - yu)\\ - C " (f+7>r ] aywuT 
" exp V (T + 7) 2 " + F+7 

If {p n \ forms the stationary distribution, this reduces to 

f *(fx + yx - yu)e~ {(+y)T + yuz 



exp 



- 1 



r + 7 

aT 1. 



awfa + yx - yu)[l - e (f+T)r ] aywuT 
(f + 7) 2 F+T 
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If, in this last expression, we let x approach 1, z approach 1, and u ap- 
proach 1, we obtain 



exp 



_ _^V6r( e — - i) _ aT 



(C) 



as the generating function E[vte~*'\ for an interval of equilibrium. 
Alternately, if instead we let x approach 1, z approach 1, and w ap- 
proach 1, we obtain (C) with u substituted for w; this implies the not- 
surprising result that for an interval of equilibrium, the two-dimen- 
sional variables {A, Z\ and {H, Z\ have the same distribution. From 
this and (C) it follows that for equilibrium (0, T), (i) .4 and H both 
have a Poisson distribution with mean aT, and (ii) the estimators L and 
hi have the same distribution. 
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